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Abstract: In order to analyze numerically inverse problems several techniques based on linear 
and nonlinear stability analysis are presented. These techniques are illustrated on the problem 
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porous media that are performed in laboratories. This is an example of the problem of estimating 
nonlinear coefficients in a system of nonlinear partial differential equations. 
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Estimation des non-linearites pour des ecoulements 
diphasiques en milieu poreux 

Resume : Afin d'analyser numeriquement des problemes inverses on presente plusieurs techniques 
basees sur l'analyse de stabilite lineaire et non-lineaire. Ces techniques sont presentees pour le 
problcmc d'estimation des mobilites et de la pression capillairc dans des deplacements dipha- 
siques unidimensionnels en milieu poreux realises en laboratoire. C'est un exemple de problcmc 
d'estimation des coefficients non-lineaires dans un systeme d'equations aux derivees partielles non- 
lincaires. 

Mots-cles : ecoulement en milieu poreux, probleme inverse, estimation des coefficients non- 
lineaires 
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1 Introduction 

Multiphase flow in porous media is modelled by a set of nonlinear partial differential equations 
equations and it provides a very good practical example for the inverse problem of estimat- 
ing nonlinear coefficients in nonlinear partial differential equations. The standard problem in 
petroleum engineering is to estimate the relative permeabilities and capillary pressure curves 
from laboratory experiments which consists of displacing a resident phase by injecting the other 
[T51 RT21 [TBI [35J [551 131 IH1 ISO] • The relative permeabilities and the capillary pressure are functions of 
the saturation of one of the phases. More recently experiments where the displacement is due to 
centrifugation were designed in order to improve the estimation of the capillary pressure function 
133 GUJ IUJ ■ Three-phase flow were also considered in [T51 HE] • In this case the relative perme- 
abilities and the capillary pressure are functions of two variables. In hydrogeology the Richards 
equation is often used and the problem of estimating its coefficients is considered in [22l [Tj . 

Without trying to give a complete review we can add to this bibliography several interesting 
contributions [T71 [331 12H and two reviews for parameter estimation in multiphase flow [TU [32] . 

In this paper we present several ingredients for a successful numerical estimation of the rela- 
tive permeabilities and the capillary pressure. In Section 2 we introduce the mathematical model 
for two-phase flow and in Section 3 we set the parameter estimation problem as a minimization 
problem. Multiscale parameterization is adressed in Section 4. Some linear analysis of the prob- 
lem is presented in Section 5 and confidence intervals are calculated in Section 6 using edgehog 
extremal solutions. Techniques for nonlinear analysis are presented in Section 7 and implemented 
numerically in Section 8. 



2 A model for a two-phase displacement in porous media 

In several laboratories core samples collected from oil fields are analyzed to determine their flow 
properties. A typical experiment consist in displacing a resident wetting fluid (subscript w), say 
water, by a nonwetting fluid (subscript nw), say oil. The displacement may be driven by injecting 
the nonwetting fluid through one extremity of the core or by centrifugal forces. These experiments 
are sketched in Fig. [T] 



Figure 1: Two-phase displacement by injection (right) or by centrifugation (left) 

Two-phase displacement is governed by a generalized Darcy's law and, in laboratory experi- 
ments, it is usually assumed to be incompressible and one-dimensional. Using the global pressure 
formulation [14j the displacement is modelled by the following nonlinear equation : 

dS dq w 

q w = -Ka(S)— + q T {t)b T (S) + q G b G (S), 

where S — S W ((Q < S < 1) is the saturation of the wetting fluid, q w its Darcy velocity, <j> is the 
porosity of the rock and K is its absolute permeability. The total flow rate qx = q w + q n w is the 
sum of the flow rate of the two phases and is given by 

DP 

q T = -Kd{S)[— - P {S)H]. (2) 
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The global pressure P [H] is given by 

P =\{p w +p nw )+l(S) 1 (3) 

with p w ,p n w the phase pressures and 7 defined below. 

H is a gravity or centrifugation function. In case of a standard displacement H{x) — gV Z(x) 
where g is the Newton constant and Z is the depth at the location x. In case of a displacement by 
centrifugation H(x) — uj 2 (r + x) where lo is the angular speed, r is the distance from the rotation 
axis to the closest extremity of the core (see Fig. [1]). The total flow rate qr is independent of x 
because of the incompressibility assumption. 

The gravity or centrifugation field qo is given by 

jr Pw + Pmv „ i A \ 

QG=K ^ H ( 4 ) 

where p w and p nw are the densities of the wetting fluid and the nonwetting fluid respectively. 

We have introduced above the coefficients a, 6t, be d, p, 7 which are functions of the saturation 
S. They relate to the relative permeability functions kr w and kr nw and to the capillary pressure 
function p c = p w — p nw through the following relations : 



kw knw t j kw , k w k nw p w Pnw 

—, Vc-, kr- ' 



Jo (6t(s) -2^' 



, 1 1 1 Pw 4~ knw Pnw 

a = K w + Knw ) P = ; — 1 7 



, kri . 
ki = , 1 = w,nw. 

where pn,i — w, nw are the viscosities of the two phases. 

The relative permeabilities kr w and kr nw and the capillary pressure p c are functions of the 
saturation which satisfy the following physical properties : 

kr w > 0, kr w increasing, kr nw > 0, kr nw decreasing, kr w {0) — kr nw (l) = 0, , . 

Pc > 0, p c decreasing. 

To equations ((2J) we add various boundary conditions depending on the experiments |36j . 



3 The parameter estimation problem 

The problem is to estimate relative permeabilities and pressure capillary functions. For this 
purpose experiments are set up so the following measurements are available : 

1. Local measurements : saturations are measured at different times tk and different loca- 
tions Xi\ 

2. Global measurements : cumulated productions Q™ — f* k (f> w (t)dt and pressure drops APk 
are measured at different times tk- 

The problem of estimating the mobility and pressure capillary functions is set as the problem of 
minimizing the error function J : 

J(kr w ,kr nw ,p c ) = SZX'^M - SZ? + ~ Qkf 

k i k 

+ 5>£(AP£-APf) 2 . (6) 

k 
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Here the superscripts m and c refer respectively to "measured" and "calculated", and w^ 1 , u>J 
and Wp are weights given to the measurements. The function J measures the difference between 
the measured quantities and that calculated with the model using the parameters kr Wl kr nw ,p c . 

The choice of parameterization of kr w , kr nw , p c will be discussed in the next section and is 
crucial for a successful estimation. If computational costs are taken into account, the choice of 
the parameterization determines also the choice of the minimization method. If the choice of 
the parameterization gives a small number of parameters, say smaller than 15, then Levenberg- 
Marquart or trust region methods [29j can be used and be cost efficient. However, when the number 
of parameters becomes large, then these optimization methods become too expensive and quasi- 
Newton methods [2] using the gradient of J calculated with the adjoint method [7[ [35] becomes 
the right choice. The drawback of this method is the difficulty to calculate the gradient and to 
implement this calculation. Therefore techniques of automatic differentiation were developped 
[25], [T9] and sophisticated software like Tapenade and Adifor are now available. 

Now we introduce some notations. We denote by A the set of admissible parameters p — 
p c ), that is parameters that satisfy properties (0. A is a subset of a set U. The direct 
mapping is the mapping ip 

(p : A dlA i — ► O (7) 

p = {kr w ,kr nw ,p c ) i — ► ip(p) = (p(u(p)) = (S c , Q c , AP C ) 

which maps A into the Hilbert space O of observations. To solve the direct problem is to calculate 
(pip) for a given parameter p. This includes solving equations ([T]), {3J with appropriate initial and 
boundary conditions. 

Let measurements z — (Q m , S m , AP m ) be given and let us write the error function J defined 
in ^ in compact form as 

Jip)=\\ <p(p)-zf w , (8) 

where |j .\\ w is the weighted norm for O with W the diagonal matrix of the weights given to the 
various measurements. The inverse problem is to minimize J over the set of admissible parameters: 

Find p S A such that J(p) = min J(p). (9) 

Of course in real life J does not vanish at the minimum because of errors in the model and in the 
measurements. 

Several questions should be adressed. Is the minimum unique ? Are there local minima ? Is 
p very sensitive to uncertainties in the measurements z ? In the next sections we shall present 
several tools that are useful to answer these questions. They are based on linear analysis (Sections 
E]and[6|) and nonlinear analysis (Sections [Jj and [8|) of stability. 

But, before, let us consider the question of parameterization, always crucial in parameter 
estimation. 



4 Parameterization 

A very common choice for parameterization of the relative permeabilities and capillary pressure 
is to use an analytical representation: 

kwi^) = 0,1V , k nw iS) — flniii(l S)^ nw , 
dP 

-^ = co + Cl S + c 2 S 2 , P C (1.) = 0. 
db 

The set of parameters to estimate is p — (a w , b w , a nw , b nw , Co, c±, C2), a set of 7 parameters. With 
such a choice it is usually not difficult to find a minimum to the minimization problem ©,(0. 

However there are many cases for which such a representation is not suitable and the relative 
permeability and capillary curves do not have such analytical shapes. To do without such a 
priori shapes a possibility is to use a discrete representation of kr w ,kr nw ,P c : p — ik rw j,j = 
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l,k 



mwj ; 



Pcj,j = l,. 



,) where k rw j,j — 1 , k rnw j , P C j are the values of kr w ,kr nw ,P c at a set 



of n s discretization points of the saturation interval (0,1). Between these points the functions 
are linearly interpolated. Note that if one uses n s — 10 equidistant saturation points, which is 
reasonable to capture nonstandard shapes, then the number of parameters is already 30. 

As reported in [3] for a standard displacement experiment, it was possible to estimate the rela- 
tive permeability curves using cumulated production, pressure drop and saturation measurements. 
However, without saturation measurement, the calculated optimal parameters were depending on 
the initial guess of the minimization algorithm. But when trying to estimate simultaneously rel- 
ative permeability and capillary pressure curves, the minimization algorithm would usually get 
stuck in a local minimum with no practical interest. 

This is the reason for the introduction of multiscale parameterization. The main idea is to 
adapt the parameterization as the minimization advances, starting with few parameters in order 
to estimate the main features of the functions, and increasing slowly their number to estimate 
their refined features. A simple example of a multiscale basis to expand a function is the Haar 
basis as represented in Fig. [H 
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Figure 2: The Haar basis 

In practice one endpoint of each of the relative permeability and capillary curves is known, 
so the continuous piecewise linear representation of the functions are uniquely defined by their 
piecewise constant derivatives which will be parameterized with the Haar basis : 

^1 = (c r w) *°(S) + EE(^)^(5), 

f=l j=l 

and similarly for — rnw ( — _ ^ — £^ — ) ]\j ^ e that with such a parameterization enforcing conditions 
dS dS 

© is simple. 

Then a multiscale optimization proceeds as follows: 

1. Minimize at scale i = 0. 

2. Augment i by 1. 
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3. Minimize at scale i using as starting point functions obtained by linearly interpolating that 
estimated with scale i — 1 . 

4. If i > i max or if scales i — 1 and i give the same estimated function, stop. 
If not go to 2. 

With such a procedure, the simultaneous estimation of relative permeability and capilary curves 
was successful while it failed when using a standard discrete parameterization 3j. Fig.[3]shows the 
progression of the multiscale parameterization for relative permeabilities and capillary pressure. 

An important advantage of multiscale optimization is that it is adaptive. It does not require 
to fix a priori the number of parameters which increases as the number of steps of multiscale 
optimization increases until the difference between two scales is not anymore significant. Therefore 
the multiscale optimization procedure determines itself an almost optimal number of parameters 
necessary to interpret the available data. 

5 Stability analysis based on the Hessian 

In this section we show how to obtain some information concerning the inverse problem using the 
Hessian of J. We linearize ip around a given parameter pq. Let us perturbate po into p with a 
small perturbation dp = p — po. Taylor's expansion gives 

ip(p) « (pfo) +tp'(po)8p. 

If we introduce Sz = z—ip(po), the observation error function can actually be written as a quadratic 
function of dp : 

J(Sp)=\\ l p'{p )6p-Sz\\ 2 w . (10) 
Then the inverse problem is to minimize this function and the minimum Sp satisfies 

H8p = ip'(p ) t W5z (11) 

with H = ip'ipoYWip'tpo) the Hessian. We note that the matrix tp'ipo), which is the jacobian 
matrix of ip at po, can be viewed as the sensitivity matrix for ip. 

Therefore solving the linearized inverse problem reduces to solving the linear system (jlip and 
this explains the importance of studying the Hessian H . 

To illustrate this we shall consider three problems of parameter estimation for which, from 
realistic data obtained in centrifugal experiments [6], we generated the measured observations 
with our simulation code. These problems are : 

1. estimating relative permeabilities while measuring saturation profiles, 

2. estimating relative permeabilities while measuring productions, 

3. estimating relative permeabilities and capillary pressure while measuring saturation profiles. 

The aim of this analysis is to study the importance of the choice of measurements for estimating 
kv w , kv nw , p c . 

Numerical results for these three problems are presented respectively in figures 01 and [51 
On each figure the parameters (before and after optimization and exact) and the observations 
(measured and calculated after optimization) are shown. The relatives permeabilities and the 
capillary pressure were discretized using the multiscale parameterization discussed in Section [4] 
with thirty parameters for each function. The saturation were observed at 6 different times in 15 
locations to give the saturation profiles which are presented. The production of the displaced fluid 
(the nonwetting fluid) was observed at the same 6 different times as the saturation profiles. 

In Figures 01 [H and E] the Hessian is also represented as a function of two variables which are 
the indices of the parameters to be estimated. They are ordered so that corresponding to the 
mobility of the wetting fluid, that corresponding to the mobility of the nonwetting fluid and that 
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Figure 3: Relative permeability and capillary pressure curves estimated with multiscale parame- 
terization 
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corresponding to the capillary pressure are in this order. For the two first experiments the singular 
values of the Hessian are also shown. Finally we drew the sensitivity of the obervations to the 
parameters that we calculated as the norms of the vector colummns of the Jacobian (ft(po), each 
colummn corresponding to the derivative of ip with respect to a parameter. 

We first observe that the parameters corresponding to small values of the saturation (smaller 
than 0.4) are not well estimated. This is not surprising since during the simulation these saturation 
are not reached (see saturation profiles). Therefore the error function is not sensitive to these 
parameters and the inverse problems is ill-posed. This is confirmed by the shape of the Hessian 
which have many coefficients which are very small and by its many zero singular values. Actually 
the direct mapping ip itself is not sensitive to these parameters as the sensitivity of the observations 
to these coefficients is zero. 

However we may notice differences between Problems 1 and 2 (observing saturation profiles 
versus observing productions). The parameter estimation works better for Problem 1 : the es- 
timated mobilities are closer to the exact ones. Another way to look at this is to compare the 
singular values of the Hessian for the two problems. We see that Problem 1 has fewer zero singular 
values of the Hessian so it is better conditionned. 

Still comparing Problems 1 and 2 we observe that the Hessian in Problem 1 is more concentrated 
along the diagonal. This indicates that for this problem the parameters are less coupled which is 
a definitive advantage when minimizing. 

These comparisons between Problems 1 and 2 give an answer to a practical question. It is 
more complicated to measure saturation profiles than it is to measure production, so are these 
efforts useful ? By analyzing the Hessian the answer is yes, and this is confirmed by numerical 
experiments SI E] - 

Finally, considering Problem 3, we observe that the saturation profiles is more sensitive to the 
capillary pressure than to mobilities. This is a confirmation of the intuition of engineers which 
designed these centrifugation experiments in order to improve the estimation of capillary pressure. 

Observation of Hessian H 

• The residual is not sensitive to parameters corresponding to small coefficients of H . 

• If H close to diagonal form, parameters are uncoupled and optimization is easier. 

• The more singular values are nonzeros, the better the conditionnement of the optimization 
problem is. 

6 Calculation of confidence intervals using edgehog extremal 
solutions 

There is a large litterature on the calculation of confidence intervals when estimating parameters 
using various methods, deterministic or probabilistic. Here, as an example, we present the method 
of edgehog extremal solutions [5H [23J 123 • 

The edgehog extremal solutions are those which correspond to parameters po + 8p satisfying 

J = o-\, || dp || 2 < erf or J <al, || Sp\\ 2 = a\ 

for given o\,o~\. a\ is the admissible maximum residual which corresponds to the error in mea- 
surements, and o~\ is the admissible maximum perturbation of the parameter. 

Assume that the matrix of weights is of the form W — wl with met and / the identity matrix, 
and introduce the residual of the linearized problem r = ip'(po)Sp — Sz. Then from equation (|10|) 
the error function J can be written as J(5p) — wr*r. 

We introduce also the SVD decomposition of the sensitivity matrix A — (p'(po) = USV 1 . 
Notice that, when the matrix W = vol, the SVD decomposition of A is closely related to the 
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Figure 4: Problem 1 : Estimating relative permeabilities while measuring saturation profiles 
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Figure 5: Problem 2 : Estimating relative permeabilities while measuring productions 
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Figure 6: Problem 3 : Estimating relative permeabilities and capillary pressure while measuring 
saturation profiles 
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spectral decomposition of the Hessian H : 

H = A l WA = VSWWUSV 1 = V diag (wsf )V* = VKV 1 

where the Sj's are the singular values of A and A = diag (A*) = diag (wsf). 

Consider a perturbation dp parallel to the eigenvector Vfc, Sp = OkYk and let us find the 
conditions it must satisfy in order to satisfy the edgehog conditions. 

Then J = wr l r can be rewritten as J = s\wa\ — 2skw(3kOk + J where (3k is the fcth component 
of the vector Sz = z — <p(po) and J =\\ Sz ||^=|| f{po) — z ||^. Therefore the edgehog condition 
J = o\ reduces to 

s\wa\ — 2skwf3kOk + J — o\ = 0. (12) 



Pk V ^ fc + CT i ~ J 

Solving for ak we obtain ak = — ± 1= ■ Since po is the calculated solution to the 

Sk VA fe 

minimization problem, (3k is small and for Sk sufficiently large - which means that the sensitivity 



to the fcth parameter is not too small - we obtain ak ~ ±— — j= — . Therefore we obtain the 



edgehog solution 



j - J 

Pe = Po + 5p « po ± - — t= — Vfc. 



/A 



This implies that, as expected, the larger Xk is, the smaller the uncertainty 5p along Vk is. 

When Xk is small, the uncertainty Sp becomes large and the edgehog condition || dp || 2 = a 



acts so 5p — ±C2 Vfc. It remains to check that J < crj. Indeed it is easy to check that in this case 
Ofc lies between the roots of the trinomial in the righthand side of equation (|12[) . 
Therefore the edgehog extremal solution associated to (Afc, Vfc) can be written as 



I- J 

PE=Po + akVk, a fe = ±min(- — = — ,cr 2 ). 

When H is diagonal a variation of the parameter Sp in the direction of Vk corresponds to a 
variation of the fcth parameter. When H is not diagonal, then the matrix H can be replaced by 
the diagonal matrices diag(/fe) or diag(^- \hij\). 

In a numerical experiment taken from [36] and whose results are shown in Fig. [3 H was 
replaced by diag(^V \hij\). The calculated minimum was J = 5.27 x 1CT 6 since we were looking 
at the synthetic example. The data for calculating the edgehog extremal solutions were eri 2 = 
2.25 x 10~ 3 which corresponds to an error of 0.005 on each saturation measurements, and o~i — 2.4 
which corresponds to a bound of 0.3 on each parameter. 

Again one can observe that for saturations smaller than 0.4 the confidence is small which is 
normal since during the experiment under study these values of the saturation are not reached. 
Actually we can observe that for values of the saturation smaller than 0.33 the extremal solution 
is determined by the edgehog condition 02 = 2.4. 



7 A geometric approach to nonlinear stability analysis 

In this section we recall results on nonlinear analysis for the problem of global minimization that 
were obtained by a geometric approach in several papers [T0l [8] . Actually in the following we 
will be in the case of a bounded admissible parameter set A in a finite dimensional parameter 
space U. We will give sufficient conditions for the problem to have a unique global solution without 
local minima and give a stability result for the global solution. These results will be applied in 
the next section to the problem of estimating the relative permeabilities. 

The following definition is devised to ensure both well-posedness and optimizability ???? of a 
nonlinear least square minimization problem. 
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Figure 7: Estimated relative permeabilities with confidence intervals calculated with edgehog 
extremal solutions 

Definition 1 The nonlinear least square minimization problem is said to be Q-well posed 

if there exists V , a neighborhood of p(A), and d, a distance on A such that, for all z G V 

1. there exists a unique global minimum p, 

2. there is no local parasitic minima for the problem (0|). fffj). 

3. the mapping z ► p is Lipshitz continuous from V C O to A. 

In order to construct such a neighborhood V we introduce some notations. For any pair po,pi 
we associate a path IT in if (A) joining the points <p(pa) and ip(pi) which is the the image by tp of 
a straight path joining pq and p\ in A. We suppose that II is twice differentiable with respect to 
its arclength s and we denote by TI'(s) and II"(s) the velocity and the curvature vectors at IT(s) 
(see Figure IH]). 

The length of the path II defines a pseudo-distance on A between any two admissible parameters 
Po and p\. We denote by 5(po,pi) this pseudo-distance on A. 

Concerning curvature we introduce not only the usual radius of curvature 



but also the global radius of curvature which is defined as follows j8|. Introduce N and N' the 
two affme subspaces normal to II at the points n(s) and n(s'), then the global radius of curvature 
Pg{s, s') between the two points n(s) and II(s') is 



This quantity is not local since it depends not only on the curve at II(s) but also at n(s'). 
Expressions for the global radius of curvature are given in Figure O One should note that 



We introduce now the set of maximum paths V = {II | po,pi S dA}, and we consider the 
worst case over one extremal path IT G V, 



PG (s,s') = d(U(s),NnN'). 



Pg(s,s') ^ p G (s',s), limpets') = p(s). 



i?(IT) = infp(s), Rg(IL) = inf p G (s,s'), 




(13) 
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Figure 8: The direct mapping tp and a maximum path. 



and for all maximum paths, 



R = inf R, R G = inf R G . 
iiev iiev 



(14) 



These numbers clearly satisfy Rc(H) < R(H),Rg < R- 



Theorem 1 Let A be bounded and Li be finite dimensional. If < Rg, then the projection on 
if (A) is Q-well posed in the neighborhood V — {z | d(z, ip(A)) < Rg} of <p{A) for the arc-length 
distance 8(po,px) on <p(A). Furthermore the following estimates hold. 
If two measurements Zq and Z\ are close enough so there exists a number d satisfying 



|| zq — z\ || +max d(zj,(p(A)) < d < Rg, 

3=0,1 



(15) 



then the following stability estimate holds for the corresponding parameters Po,Pi obtained by 
solving the associated least square problems : 



x( x ^ g(g) n || 
d(Po,Pl) < T^rprr \\ Z - Z\ || 



R(U) - d 

where II is the path connecting f(po) and ip(pi). 



(16) 



Inequality (fT!)|) is a stability result for the arc length distance f{A). Depending on the hypothesis 
made on <pf(p), it will imply two stability estimates on || po — pi \\ given below, the second one 
being sharper than the first one. 



Q-well posedness: 

Assume that 



There exists a m > Osuch thata m || q \\<\\ f'{p)q || for allp £ A, q G U. 



(17) 



Then we have 



Combining this inequality with estimate (|16J) we obtain 

1 R 



Hpo,Pi) = / Wipa +t(p -pi))(po ~Pi)\\dt > a m \\ p a -pi || . 
Jo 



Po-Pi \\< 



R-d 



II z o ~ zi || 



(18) 



and problem ©,((8]) is Q-well posed. 
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Figure 9: The global radius of curvature 
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Directional stability: 

To improve the above estimate, we give now a directional estimate around a specific point, say 
Po- We write p\ E A as p\ = po + hv with h —\\ po — p\ || and v a unit vector. We now assume, 
instead of (fT7]) : 

There exists a m (po 1 v) > such that a m (po,v) <|| ¥>'(po +tv)v ||, for allpo + tv E A. (19) 
This is a less demanding condition than (|17p since a m < a m (po,f). Then we have 

<5(Po,Pi) = / \\<p'(po + thv)hv\\dt > a m (p ,v)\h\. 



Combining with (|16p we obtain now 

II Po ~ Pi ||= N < r fl(po 'P l} || z - Zl || (20) 

a m {pa,Pi) R(Po,Pi) - a 

where R(po,pi) = R(H) is just a notation stressing the dependence of the smallest radius of cur- 
vature along II between po and Pi ■ 



8 Implementation of the geometric nonlinear analysis 

In this section we show how to use Theorem [T] in practice to estimate uncertainties in the param- 
eters from uncertainties in the measurements. We consider a two-phase displacement where we 
estimate relative permeabilities of the form 

kr w (S) = S a , kr nw (S) = (l-S) b . 

Here a and b are the parameters to estimate. The constraints that we impose on them are 
1 <a<3,l <6<3so„4 = (1,3) x (1,3). The observations are saturation profiles. Saturations 
are measured at 6 different times in 5 different locations (30 measurements) in a first case, and in 
15 different locations (90 measurements) in a second case. We set up the experiments so that we 
know the optimal parameters : d = b = 2. 

To calculate estimates for Rq , R, a m we proceed as follows : 

1. Choose a sample of maximal paths V* in A that we assume is large enough to represent V 
the set of maximal paths. An example is given in Fig. [5] 

2. Discretize each paths with a set of points. 

3. Calculate at these points II, II', II", these derivatives being made with respect to arc length 
and being calculated for instance by finite differences. Remember that each calculation of 
II at one point requires a solution of the direct problem. 

4. Calculate for each path II, i?c(n), i?(n), a m (II). For that we use equations (fT4|) where we 
replace the infimum over all points of a path by that over the set of discretization points 
and the infimum over V by that over the subset V* of V . a m (II) is estimated by using 
finite differences with the points discretizing II. The results are given in Table [T] for 30 
measurements and for 90 measurements. From these results we obtain approximate values 
Rq, R*, a* m of i?, a m for 30 measurements, 

R* G = R G {AG) = 1.11 x 10~ 3 , iT = R{AG) = 1.11 x 10" 3 ,a r * n = a{BF) = 0.077, 

and for 90 measurements, 

R* G = R G {EG) = 1.70 x IQ- 2 ,R*R{EG) = 1.70 x 10~ 2 ,c4 = a(AG) = 0.28. 
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G F E 

3 _ , , , 




A B C 



Figure 10: A sample of maximal paths for A = (1,3) x (1,3). 



n 


Ra(n) 


R(U) 


a™ (II) 


R G (n) 


R(U) 




AC 


2.34 x 


2.34 x 10^ 


0.267 


6.44 x 10~* 


6.44 x 10~ 2 


0.545 


AD 


2.12 x 10~' z 


2.12 x W~' z 


0.244 


7.59 x W~' z 


7.59 x 1Q- Z 


0.457 


AE 


i.7i x ltr^ 


1.71 x 10^ 


0.205 


4.95 x 10^ 


4.95 x 10~ 2 


0.392 


AF 


1.39 x W~' 2 


1.39 x 10~* 


0.121 


3.56 x W~' 2 


3.57 x 10~ 2 


0.297 


AG 


1.11 x 10 a 


1.11 x 10 a 


0.520 


3.07 x W~' z 


3.07 x 10~ 2 


0.280 


BD 


2.20 x 10^ 


2.20 x 10^ 


0.236 


3.92 x 10^ 


3.92 x 10~ 2 


0.422 


BE 


9.46 x 10~ 3 


9.52 x 10~ a 


0.140 


2.70 x 10^ 


2.71 x 10~ 2 


0.360 


BF 


7.44 x 10~ 3 


7.44 x 10~ 3 


0.077 


2.88 x W~' z 


2.88 x 10~ 2 


0.348 


BG 


2.75 x 10~ 3 


2.75 x 10~ a 


0.138 


2.68 x 10^ 


2.68 x 10~ z 


0.375 


BH 


1.98 x 10^ 


1.98 x 10~* 


0.189 


6.09 x 10^ 


6.09 x 10~ 2 


0.500 


CE 


4.90 x 10 -3 


4.90 x 10~ d 


0.085 


2.92 x 10^ 


2.94 x 10~ 2 


0.414 


CF 


1.38 x 10^ 


1.39 x 10^ 


0.193 


4.67 x 10^ 


4.67 x 10~ z 


0.514 


CG 


4.01 x 10~ 3 


4.50 x 10~ a 


0.190 


2.20 x 10^ 


2.20 x 10~ 2 


0.421 


CH 


2.15 x 10~' z 


2.15 x W~' z 


0.223 


6.07 x 10^ 


6.07 x 1Q- Z 


0.520 


DF 


1.92 x 10~' z 


1.93 x W~' z 


0.266 


6.04 x W~' z 


6.05 x 10~ 2 


0.606 


DG 


5.43 x 10~ 3 


5.91 x 10~ a 


0.227 


2.04 x 10^ 


2.09 x 1Q- Z 


0.437 


DH 


1.78 x 10^ 


1.79 x W~' z 


0.236 


6.28 x 10^ 


6.29 x 1Q- Z 


0.476 


EG 


6.58 x 10~ 3 


6.58 x 10~ a 


0.244 


1.70 x lO'* 


1.70 x 10 ^ 


0.403 


EH 


2.61 x 10^ 


2.61 x 10^ 


0.202 


4.17 x 10^ 


4.18 x 10~ z 


0.365 


FH 


1.85 x 10~' z 


1.85 x 10 _;a 


0.157 


3.43 x 10^ 


3.43 x 10"^ 


0.280 




30 measurements 


90 measurements 



Table 1: Values of i?c(n), a m (II), -R(II) for all paths II e V* when using 30 and 90 measurements 
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\AS\ 


II Sz || 


d = 2 || Sz || 


II s P || 


30 measur. 


7.3 x 10" 5 


0.0004 


0.0008 


0.0191 


90 measur. 


7.3 x 10 _b 


0.0007 


0.0014 


0.0027 



Table 2: Uniform estimates on the error 5p on the parameter for a given error \AS\ on a saturation 
measurement. 

We can now apply Theorem [T] 
Q-well posedness 

In the case of 30 measurements, if the error on the measurements Sz is such that 

II Sz ||< R* G = 1.11 x 1CT 3 , 

which corresponds to a 2 x 10~ 4 error on each saturation measurement, then it follows that the 
measurement z lies in the neighborhood V. We see that (fTT)) is satisfied with a m = 0.077 so that 
the nonlinear least square problem (O,© is Q-well posed. 

Similarly, in the case of 90 measurements, if the error on the measurements Sz is such that 

|| Sz \\< R* G = 1.7 x 10~ 2 , 

which corresponds to a 1.8 x 10~ 3 error on each saturation measurement, then it follows that the 
measurement z lies in the neighborhood V. We see that (|17p is satisfied with a m = 0.28 so that 
the nonlinear least square problem ©,(11]) is Q-well posed. 

One can notice that, as expected, increasing the number of measurements allows for larger and 
larger errors on saturation measurements. Practically, even with 90 measurements, the precision 
required for the saturation measurements is difficult to achieve. 

Stability 

Given po = p the estimated parameter, which is a minimizer of the nonlinear leat square 
problem we can use the stability result of theorem [T] to estimate the uncertainty Sp in the 

following way. Denote pi — Pq + Sp and assume that we know the uncertainty |A)S| = 7.3 x 10 -5 
on one saturation measurement. This uncertainty corresponds to an uncertainty || Sz | on the 
vector z of saturation measurements (see table [5]). We take d = 2 || Sz || so that d < Rq and 
hypothesis (fT5| is satisfied. Then inequality (flgf . when replacing R and a m by R* and a* m gives 
the uniform bounds on || Sp || given in table [21 These estimates on || Sp || can be represented by 
the domains of uncertainty shown in Fig. [TTJ circles centered at po of radii || Sp \\ . Note that these 
domains of uncertainty do not actually depend on the value of the optimal parameter p$. 

If, instead of inequality ([T5]) , we use inequality (j2"0"|) to study the uncertainty with respect to 
saturation measurements, we proceed in the following way to build the domain of stability shown 
in Fig. fT2l around the calculated minimum p = (2, 2). There are 6 paths of V* going through p . 
On each of this path we look for the parameter p\ the furthest from po satisfying inequality ([20]) . 
The domain of uncertainty is now an irregular polygon whose shape depends on pq. Comparing 
the scales of Figs. [TT1and[T2l we observe that the domains of uncertainties obtained from 1st order 
directional estimate ([2Tfl) are significantly smaller than those obtained from uniform estimate ([T8]) , 
confirming that estimate PO]) is sharper than estimate ([T5]) . 

We finally remark that in any case, having more measurements reduces the domains of uncer- 
tainty. 

9 Conclusion 

When estimating the relative permeability and capillary pressure functions in two-phase displace- 
ment experiments, we showed that much information about stability and uncertainty on the esti- 
mated parameters can be obtained from the Hessian. 



RR n° 6892 



20 



J. Zhang, G. Chavent & J. Jaffre 



Figure 11: Domains of uncertainty from uniform estimate (|18p . 




1.992 1.994 1.996 1.998 2 2.002 2.004 2.006 2.008 

a 

Figure 12: Domains of uncertainty from 1st order directional estimate (|20j) . 
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Furthermore we showed also how to use in practice geometric nonlinear analysis tools to claim 
optimizability and to construct domains of uncertainty. 
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